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The development of a wavelet-based feature extraction technique specifically targeting FOD- 
event induced vibration signal changes in gas turbine engines is described. The technique performs 
wavelet analysis of accelerometer signals from specified locations on the engine and is shown to be 
robust in the presence of significant process and sensor noise. It is envisioned that the technique 
will be combined with Kalman filter thermal/health parameter estimation for FOD-event detection 
via information fusion from these (and perhaps other) sources. Due to the lack of high-frequency 
FOD-event test data in the open literature, a reduced-order turbofan structural model (ROM) was 
synthesized from a finite element model modal analysis to support the investigation. In addition to 
providing test data for algorithm development, the ROM is used to determine the optimal sensor 
location for FOD-event detection. In the presence of significant noise, precise location of the FOD 
event in time was obtained using the developed wavelet-based feature. 

1.0 Introduction 


Ingestion of foreign objects by a turbofan engine, particularly birds and ice, is one of the leading 
causes of uncontained rotor events in commercial jet transport. There are numerous cases of bird 
ingestion contributing to accidents, some of them fatal [1]. Specific procedures have been 
developed for situations where ingestion is suspected, and forensic analysis has repeatedly shown 
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that they would have been the appropriate actions to take in cases where rejected takeoffs 
motivated by bird ingestion resulted in accidents. These are the worst situations; in the vast 
majority of cases, bird ingestion does not affect the safe outcome of a flight and may, in fact, go 
unnoticed by the flight crew. However ingestion does pose a risk of Foreign Object Damage 
(FOD), and thus is important to identify, if possible. When a multi-engine aircraft flies through a 
flock of birds, potentially damaging more than one engine, the pilot needs to understand the status 
since his resulting actions may be different depending on the number of engines involved. Even in 
a case where there is no apparent damage to an engine as observed from the cockpit, latent effects 
(e.g. cracks that can be propagated by high cycle fatigue) may be present that require maintenance 
action [2], 

The critical consequence of foreign object ingestion is engine surge, potentially resulting in the 
loss of power. The flight crew can recognize foreign object ingestion through a combination of 
instrument readings and sensory cues. These include such symptoms as a thud or bang, a fire 
warning, a visible flame coming out of the engine, vibration, yaw of the airplane caused by thrust 
imbalance, high Exhaust Gas Temperature (EGT), change in the spool speeds, smoke/odor in cabin 
bleed air, and Engine Pressure Ratio (EPR) change. It is important to note that for impact-type 
FOD (due to ice, birds, runway debris, etc.), the damage is primarily to the fan and front part of the 
engine, with the extent of the damage determined by the geometry, angle of impact, hardness, 
relative speed, etc. of the object. 

One key reason for detecting FOD automatically is that most such events occur close to the 
ground when the flight crew is occupied flying the plane. Only after the crew stabilizes the aircraft 
at a safe altitude should they take action on the engine. Once the aircraft is stable, the reduced 
workload in the cockpit environment provides better circumstances under which the crew can 
analyze the situation, and they would benefit greatly by having full knowledge of the engines 
involved and the likelihood of the event based on data. 

Given that many of the potential symptoms of FOD are not unique to that type event, an 
automatic system for FOD detection should take information from multiple sources to provide 
confidence in its diagnosis. Just like a pilot does, this system must fuse information in a way that 
provides a measure of the likelihood of foreign object ingestion in order to determine any 
corrective action. Additionally, since several of the potential symptoms are described in terms of 
human senses, alternate information sources need to be developed. These sources should utilize the 
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standard engine sensor suite (or a very similar suite) to address issues such as certification and 
retrofit. Naturally processing requirements may vary from current systems simply due to the fact 
that additional on-line signal processing would be required, but the impact would be minimal since 
such a system would not be flight critical. 

For approximately two decades techniques based on the Kalman filter algorithm have been 
applied to turbofan engine diagnostics [3,4]. Specifically, analytical (or virtual) measurements 
provided by Kalman filter estimators have been used for detection of engine degradation via 
estimation of a set of health parameters, which are generally not measurable themselves (e.g., 
compressor efficiency) and are calculated via knowledge of measurable quantities. The degradation 
monitored may be gradual in nature, e.g., worn components that increase internal clearances, 
resulting in decreased component efficiency. It may also be abrupt as is the case when foreign 
objects are ingested by the engine. Changes in component efficiencies, high and low spool speeds, 
and thermodynamic parameters have been determined/observed during foreign object damage 
(FOD) events [3,4]. In some cases, these parameter changes alone may not be conclusive proof that 
a FOD event has occurred. Absent from these investigations is the incorporation of structural 
vibration signals as a means to aid in positive identification of a FOD event. 

To successfully apply available sensor fusion techniques, the same event should be detected 
using sensors that have significantly different physical characteristics (e.g., thermocouples 
compared to accelerometers) and rely on measuring physically different parameters (e.g., 
temperature and acceleration) [5,6]. It would appear that fusing health parameter estimates with 
structural response information acquired during a FOD event could provide conclusive evidence 
that a FOD event had occurred. Furthermore, fusion of virtual measurements and vibration signal 
features may enable discrimination between a normal aircraft maneuver/operating condition and a 
more severe FOD event even though both events appear to have similar looking health parameter 
data profiles (due to inherent limitations of the diagnostic model used). Similarly, vibration signals 
alone may possess FOD-like characteristics depending on the maneuver e.g., a rapid change in 
thrust demanded from the engine would result in transmission of temporary impulse-like force 
imbalances through the engine structure. The combination of traditional health parameter 
estimation and of state-of-the-art signal processing techniques applied to structural vibration 
signals (e.g, wavelet analysis [7,8]) could provide the “finger print” necessary to positively identify 
a FOD event. 
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This paper describes the development of a wavelet-based feature (WF) extraction technique 
specifically targeting FOD-event induced vibration signal changes in gas turbine engines. It is 
envisioned that the wavelet-based feature will be fused with measurements or estimates of engine 
thermal parameters, themselves having FOD-specific response characteristics, to provide a high- 
confidence FOD event diagnosis. Successful automated identification of FOD events would be 
especially useful in unmanned (autonomous) vehicle applications, where the absence of the pilot 
precludes the use of human monitoring of engine conditions. 

To the authors’ knowledge, engine-specific vibration signatures associated with actual FOD 
events e.g., rotor frequencies excited due to foreign object impact on the fan disk, is not available 
in the open literature. Potential sources of high-frequency structural data, specifically high-fidelity 
finite element models of aircraft engines, do not easily lend themselves to sensor placement and 
diagnostic techniques based on well established linear system theory. As well, these models are 
typically computation-time intensive due to the vast amounts of time-dependant spatial data they 
produce and do not lend themselves for direct application to the initial stages of control and 
diagnostic (CDS) development and testing. To facilitate WF development, a reduced order model 
is synthesized to provide the structural response signals indicative of a FOD event. The low-order 
structural model is also used to determine the optimal accelerometer locations for feature 
extraction. This requires a linear model in state-space form to facilitate the calculation of an 
Observability Grammian matrix [9] that would determine whether a mode excited by a FOD event 
could be detected by a particular sensor arrangement. 

2.0 Development of reduced-order turbofan rotor structural models 

Finite element structural analysis codes provide invaluable information concerning structural 
deformation, stress intensities and modal behavior. However, finite element codes tend to be 
cumbersome and time consuming to work with when used to design and test control and diagnostic 
system (CDS) algorithms, especially when performing multiple design iterations. State-space 
models, on the other hand, are the form of choice for many signal processing and Multiple Input 
Multiple Output (MIMO) controller applications. A method for reduced-order modeling of 
turbofan rotor structures is presented which uses output from a finite-element code in the form of 
natural frequencies and eigenvectors to create a state-space model suitable for use in CDS 


NAS A/TM— 2004-2 13118 


4 



development and testing. The resulting state-space model can easily be ported to control system 
development software such as the MATLAB/SIMULINK® suite of software tools [10]. The 
following presents the reduced-order model development process as applied to the modeling of a 
generic turbofan rotor, ultimately used for generating a simulated FOD event structural response. 

For the purpose of developing reduced-order structural models for application to CDS 
development, a dynamic model based on the dominant modes of vibration of a turbofan rotor- 
bearing system is considered adequate. The modal behavior characterized by the model should be 
representative of that presented in the literature for turbofan engines [11,12]. For the present study 
a specialized rotor dynamics finite element code, DyRoBes [13], is used to simulate the structural 
response resulting from unbalance and foreign object damage. This code is suitable for the present 
study since it has the capability to model critical features of an engine system and can sufficiently 
simulate engine dynamic response, particularly the response at sensor locations. Being a special- 
purpose Finite Element code i.e., designed specifically for rotor dynamic analysis and users with a 
basic knowledge of FE methods, it is also relatively straightforward to use. 

The structural model used for this study was derived from a generic engine model created to 
develop an industry standard for engine-airframe structural analysis. The model is of a dual spool 
rotor mounted on flexible supports (Figure 2.1). The flexible supports are characterized to simulate 
the static engine structure and wing. 

The first spool is representative of the fan and low pressure turbine stages while the second 
spool represents the high pressure compressor and turbine. An interstage bearing connects the two 
spools together. The low pressure spool is 112 inches long and consists of 18 shaft finite elements. 
A large disk is located at the left end to represent the fan and several smaller disks are located on 
the right end for the low pressure turbine. The entire shaft weighs 2494 pounds. The high pressure 
spool is 80 inches long and consists of 16 shaft finite elements. There are eight disks on the shaft’s 
left end representative of a compressor. One disk is located on the shaft’s right end representative 
of the high pressure turbine. The high pressure shaft weighs 1280 pounds. Each shaft is supported 
on bearings as shown in the figure. The bearings in turn are supported on flexible supports 
representative of the support flexibility provides by the static engine structure. The bearing and 
support stiffnesses are assumed to be 1.0 x 10 lb/in and 1.0 x 10 Lb/in respectively. 

The transient response from rotor imbalance and FOD of the turbofan DyRoBes model is 
shown in Figure 2.2. A typical rotor vibration sensor signal measured during a FOD event (with a 
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sensor located at one or more of the bearing locations shown in Figure 2.1) would consist of the 
response of a dominant low-frequency lateral mode excited from an impulsive moment at the fan 
disk / rotor shaft interface (about the X-axis in Figure 2.3) due to the initial impact of the foreign 
object upon the fan disk. Higher frequency lateral modes would also be excited due to the high 
frequency content of the foreign object impact and would appear in the sensor 



signal immediately after the event. Possible longer lasting effects due to resulting permanent blade 
damage i.e., an imbalance force occurring at a frequency corresponding to once per revolution of 
the rotor, may also appear. Rotating equipment tends to have a relatively small degree of inherent 
imbalance after manufacturing [11], which will also be included in the reduced order structural 
model for a “steady state” response. To the authors’ knowledge, high-frequency accelerometer 
data from actual turbofan engine FOD events does not appear in the open literature, thus the 
hypothesized sequence of events are predominantly based on the authors’ experience and intuition. 
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The following section describes the method used to develop a Reduced Order Model (ROM) that 
closely approximates the actual rotor response but is more amenable to real-time assessment of 
feature extraction techniques used for identifying FOD. 



Figure 2.2. Example FOD event FEM model vibration response characteristics. Calculated 
displacements are functions of the imbalance magnitude and foreign object impact 
characteristics (e.g., relative speed, size, etc...). 


2.1 Dynamic Representation of Rotor Bearing Systems in the State Space 


Dynamic systems may be represented as transfer functions or in state-space [14]. For CDS 
development the state-space form is often preferred and will be the form used for the reduced order 
model synthesized from DyRoBes FEM modal analysis data. For a multiple degree of freedom 
(MDOF) translational system the equation of motion in physical coordinates is: 


[M|q}+[cfq}+[Klq} = {F} 2.1 
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with the damping matrix C typically assumed to be proportional to the mass and stiffness matrices 


[C] = «[M] + P[K] 


2.2 


The dynamics of rotating systems are complicated by the presence of gyroscopic effects 
[1 1,15]. Expanding equation 2.1 in terms of translational and rotational degrees of freedom and 
including gyroscopics yields: 
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2.3 


Unlike a non-rotating system, the damping coefficient matrix in Equation 2.3 is a function of 
the shaft rotational speed, Q i.e., gyroscopic effects result in the rotor’s lateral natural frequencies 
and mode shapes being a function Q . To determine the natural frequencies and mode shapes, 
assume that the free response of the system is sinusoidal in nature i.e., 

r f , iCD t 

{q(t)} = {<p}e n 2.4 

where {cp}is the mode shape (eigenvector) corresponding to co n , the natural frequency (eigenvalue). 
N is the number of modes included in the model. Inserting Equation 2.4 into the unforced, 
undamped version of Equation 2.3 gives 

+ K0 n [G]+[K])[<p}=0 2.5 
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where the gyroscopic matrix G is given by 


0 0 0 0 

0 0 0 0 

0 0 0 -2I a Q 

0 0 2I 0 Q 0 


2.6 


A trivial solution to Equation 2.5 is obtained with{cp}=0. For a nontrivial solution 


(-co n 2 [M] + /Q) n [G]+[K])=0 


2.7 


must be satisfied. To obtain a reduced order linear model with coefficients that are independent 
of Q , i.e., with constant natural frequencies co ni (i=l..N) and modal damping, it is assumed that 
the natural frequencies and modes shapes calculated via Equation 2.7 at the rotor critical speeds 
adequately account for the effects of gyroscopics over the operating range of interest. For the rotor 
described in Section 2.0, the co ni and cp ; are calculated using the DyRoBes FEM code. The mode 
shapes cpj obtained are orthogonal with respect to the mass and stiffness matrices and provide the 
following normalization 


fe} T [Mfo,}=l 2.8 

{<P i } T [Kl<P i }=^ i 2.9 

Given the assumptions described above, the equations of motion for the linearized rotational 
system take the form of Equation 2.1. For lightly damped systems the matrix C is typically 
assumed to be diagonal [16]. With the eigenvectors mass-normalized and [®] = [cpj •••cp N ], 
the coordinate transformation 


{q(t)}= [®]{r|(t)} 


2.10 
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is inserted into Equation 2.1 and pre-multiplied by ® T . The dynamics of the system in modal 
coordinates are 


{ii}+ diag[2^co ni ]{q}+ diag[co^ ]{r)}= ® t {f} 


2.11 


with the corresponding state-space form 


-diag 
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2.12 


and the modal damping ratios defined as [16] 

5,= C 


c 


i critical 


2.13 


For the present study, the coordinates of interest (Figure 2.3) are 

w-B 

at specified points (i.e., node locations in Figure 2.1) in the rotor-bearing system. 

The above presents a simplified formulation for the equations of motion of a rotational system 
which, in a linearized sense, include the gyroscopic effects. These effects are accounted for in 
DyRoBes’ analyses. Modal damping ratio estimates are based on whirl analysis log-decrement 
calculations at the rotor-bearing system critical speeds and the nominal speed of rotation. 
Comparison of the reduced-order structural model frequency response (to an assumed imbalance) 
and transient response to the FEM-generated data will provide a measure of the fidelity of the 
reduced-order model to appropriately describe a FOD-event induced structural response. 
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Figure 2.3. Coordinate system used in the FEM and reduced order state space structural models. 
Cyclic imbalance forces are perpendicular to the Z-axis. Impulse moments due to a FOD event are 
about the X-axis. 

2.2 Foreign Object Damage Events modeled using Reduced Order Rotor Models 

from FEM output data. 

The state space model presented in Equation 2.12 requires the modal or natural frequencies of 
the system and the eigenvectors (i.e., modal displacements) at specified locations on the rotor or 
bearings / supports. For the purposes of FOD detection algorithm development, the locations of 
interest would be those corresponding to possible vibration sensor locations i.e., nodes 
corresponding to the bearings and supports as shown in Figure 2.1. Commercial FEM codes can 
readily provide this information [13]. Transformation from modal coordinates to physical 
coordinates is performed via Equation 2.10, with the matrix ® formed by columns corresponding 
to the eigenvectors for each frequency being considered. To transform external forces on selected 
locations on the shaft (or supports) from physical to modal coordinates, the external force vector is 
pre-multiplied by ® T . 

Physical (continuous) systems are inherently infinite order. Thus FEM, being lower order than 
the physical system being modeled, are implicitly inaccurate with FEM being more accurate in the 
low frequency range. The frequencies of interest are from 0 to 20,000 rpm, with the nominal 
rotational speed of the low pressure shaft being 9000 rpm. This frequency range will provide the 
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required response due to fan imbalance, as well as the impulse response characteristics due to FOD 
impact upon the fan disk. For adequate imbalance response over the full range of operating speeds, 
the natural frequencies and mode shapes obtained from the FEM correspond to the so-called critical 
speeds of the shaft [1 1], To attain proper response to a FOD event, several of the lateral modes 
calculated via an FEM whirl analysis at a specified rotor speed are included in the model as well. 
The number of lateral modes included is determined by the performance of the ROM with respect 
to the FEM. Once again, the goal of the present study is to develop structural models to test real- 
time diagnostic algorithms, with these models providing acceptable response characteristics for the 
forcing functions of interest (imbalance and foreign object impact) while maintaining the model 
order low. The reduced-order model development technique described above can be applied to 
FEM developed using codes intended for more elaborate analyses. 

The reduced order dynamic model was implemented in MATLAB®, using the companion 
software package SIMULINK®, in state-space form. Figure 2.4 depicts the SIMULINK® 
structural model. Gas turbine rotors tend to have an inherently small degree of imbalance, which is 
modeled as an input corresponding to a variable amplitude sine wave with a frequency equal to the 
high and / or low speed rotor speeds. Imbalance force amplitude is frequency dependent according 
to the following relation 


^imbalance — G 


2.15 


where m is the total mass of the rotating component, e is the effective eccentricity in inches, and Q. 
is the shaft rotational speed in rad-s _1 . If a FOD event results in additional imbalance, the 
eccentricity is increased by an additional amount. For the present study, the fan is assumed to have 
an eccentricity (i.e., inherent imbalance) of 0.001 inch. 

A FOD event is modeled as an object hitting a point along the radius of the fan, resulting in a 
moment input at the node on the rotor corresponding to the fan location (Node 1 in Figure 2.1). The 
magnitude of the input moment is determined as follows. Given a foreign object of mass m FO with 
an initial linear velocity relative to the fan disk v FO along the z-axis (Figure 2.3), conservation of 
angular momentum may be applied to determine an estimate of the initial angular velocity 0 of the 
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Figure 2.4. MATLAB/SIMULINK implementation of the reduced order turbofan 
rotor/bearing system structural model. 
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fan disk about an axis perpendicular to the axis of rotation of the rotor (X-axis in Figure 2.3) due to 
the impact i.e.. 



where I e is the moment of inertia of the fan disk about the X axis (Figure 2.3) and r is the distance 

of the impact from the axis of rotation (Z axis in Figure 2.3). The corresponding moment applied to 
the rotor at node 1 from the impulse load on the fan disk is determined by the rate of change of 
angular momentum, given a specified time interval for the impact to occur At . 


F* 



2.17 


The above treatment provides a straight-forward means for determining the impulse-load related 
torque input to the rotor during the initial stages of the FOD event, and is considered to be adequate 
within the level of complexity of the model needed for real-time diagnostic system development. 

As mentioned previously, information regarding specific FOD events does not appear in the open 
literature. The FOD event impact time interval, At , is dependant upon random factors such as 
foreign object geometry, location and angle of impact, foreign object hardness, relative speed, etc. 
and is estimated based on the authors’ “intuitive feel” for an event that would cause damage. 

Figure 2.5 presents a frequency response comparison of 4, 10, and 25 degree of freedom 
(DOF) reduced order models with the frequency response of the 154 DOF DyRoBes FEM (refer 
to Figure 2.1 for node locations). Figures 2.6 and 2.7 present the displacement response of the 
10 DOF ROM and FEM to a combined FOD event - imbalance input. Figure 2.6 presents the 
response at bearing / support location 37 to a “large” FOD event (m F o = 2 lbm) and a “small” FOD 
event (m F0 = 0.25 lbm). Figure 2.7 presents the response at bearing / support location 39 to a 
“large” FOD event (m F o = 2 lbm) and a “small” FOD event (m F o = 0.25 lbm). The imbalance input 
is calculated using Equation 2.15 with an eccentricity of 0.001 inches concentrated at the fan node 
(node 1 in Figure 2.1). The magnitude of the FOD event moment input at node 1 is calculated using 
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Equation 2.16 with a radius of 20 inches, specified mass, a At of 0. 1 seconds, and a relative 
velocity of 300 mph. To provide adequate FOD event response characteristics, the first lateral 
mode from an FEM whirl analysis at 9000 rpm was included. The 10 DOF model provides 
acceptable correspondence to the FEM and is considered sufficient for the present study. 
Differences between the reduced order model and the FEM are due to the assumptions made, i.e., 




Figure 2.5. Steady state imbalance response comparison of 4, 10, and 25 degree 
of freedom (DOF) reduced order models for two different bearing/support 
locations. Imbalance has an eccentricity of 0.001 inch located at the fan. 
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Figure 2.6. Response at bearing/support location 37 to (a) a “large” FOD event 
(m F0 = 2 lbm), and (b) a “small” FOD event (m F0 = 0.25 lbm). 
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Figure 2.7. Response at bearing / support location 39 to (a) a “large” FOD event 
(m F0 = 2 lbm) and (b) a “small” FOD event (m FO = 0.25 lbm) (b). 


NAS A/TM— 2004-2 13118 


17 


the ROM uses constant modal damping ratios and natural frequencies over the full range of 
operating speeds. 

For both locations, the change in the vibration signal due to the large FOD event is easily 
discernable from the nominal vibration signal due to inherent imbalance alone. For the small FOD 
events, there is less of a difference from the normal operating characteristic. For smaller FOD 
events with the presence of sensor or process noise (which would be the case in an actual engine) 
the event could easily be masked. Thus, the signal analysis technique employed to detect such an 
event must be capable of identifying the characteristic vibration signature in the presence of 
process and sensor noise. 

3.0 Optimal Sensor Placement for FOD event detection. 

Signal processing techniques that provide information concerning the internal modes (i.e., 
states ) of a system require that these state variables be observable using the available sensor suite 
[17]. For FOD-event detection, the vibration sensor (e.g., accelerometer) signal must contain the 
modal response of the rotor in the frequency range of interest. The modes of a system are 
observable if and only if the matrix 


0 = 


C 

CA 


C A 11-1 


3.1 


known as the observabilty matrix, has rank n, where n is the dimension of the state space of the 
linear time invariant system under consideration. For the 10 mode ROM representing the rotor- 
bearing system shown in Figure 2.1 (with nodes 37 through 41 as candidate accelerometer 
locations) the system is completely observable. To determine the optimal (or worst case) 
accelerometer location for FOD event detection, a measure of the degree of observability of a 
system for a given sensor arrangement must be determined. One method suggested in [17] requires 
calculation of the Observability Grammian to determine which states, starting from a specified set 
of initial conditions, have an influence on the output energy E(x(0)) of the system over all time. 
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Consider the following measure of the output energy at some initial condition x(0) for the system 
of Equation 2.3 under observation with no external input 


E(x(0)) = J||y(t)|| dt = J y(t) T y(t)dt = x(0) T Q x(0) 3.2 

0 0 

where the Observability Grammian is defined as 

Q = Je A t C T Ce At dt 3.3 

o 

For a stable system, the existence of Equation 3.3 is guaranteed and Q may be found by solving 
the following Lyapunov equation 


A t Q + QA + C t C = 0 3.4 

For a reduced number of accelerometers (reduced relative to the initial set i.e., accelerometers 
37 through 41 in Figure 2.1) the system is observable if and only if Q is positive definite. With 
observability established for the reduced sensor suite system, Equation 3.2 is used to determine the 
output energy contribution of each state variable. Large output energies correspond to easily 
observable state variables (modes). The optimal suite of sensors would naturally be that with the 
(relatively) largest output energies for the modes needed to identify a FOD event. Consider the 
calculation of Equation 3.2 for locations 37 and 40. With the entry in the initial state vector (i.e., 
x(0)) for a specified mode set to 1.0 and the initial conditions for each of the other modes set to 0.0, 
E(x(0)) for each mode observed separately is simply the diagonal entry of Q corresponding to each 
mode. Thus, a judgment as to the optimal placement of the accelerometers may be made based on 
the diagonal entries of Q alone. For locations 37 and 40 
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diag[Q] 

Node 37 — 


4. 121 le- 003 
5.1742e-003 
2.3575e-005 
9.3601e-007 
4.4442e - 006 
3.9526e-004 
2.51 lie -003 
3.6104e-004 
5.0943e-003 
9.8736e-004 


3.6160e - 007 
3.1261e - 005 
1.4738e-005 
6.9672e - 003 
4.4159e-004 
4.3813e-004 
1.6479e-004 
2.9526e - 004 
2.2466e - 008 
1.0822e - 004 


diag[Q] 

Node 40 


3.5 


which, due to the higher observation energies associated with the first three modes, implies that 
placement of an accelerometer at location 37 may be more appropriate for detecting a FOD event. 
Upon performing the calculation of Equation 3.2 for nodes 37 through 41 as candidate 
accelerometer locations, node 37 was determined to be the preferred location with node 40 being 
the least preferred. Thus, the analysis (i.e., feature extraction) technique selected should be applied 
to the signal generated by “accelerometer 37” in order to have the highest likelihood of detecting a 
FOD event. However, placement of a sensor at the optimal location may not be physically 
attainable, thus it is desirable that this same technique should be robust enough to detect the event 
at a suboptimal (possibly the most limiting) location. For some systems this conclusion may be 
drawn directly from physical intuition. For more complex systems (which may not easily lend 
themselves to physical intuition) the analytical treatment presented above provides a quantitative 
measure for optimal sensor placement. 


4.0 FOD Event Detection Via Rotor Bearing Vibration Signal Analysis 

Consider the simulated bearing accelerometer signals shown in Figure 4.1. To demonstrate the 
robustness of the WF extraction technique, the analyses described below were performed on the 
signal from the most limiting accelerometer location on the rotor bearing system considered, which 
was determined from the sensor placement analysis of section 3.0 to be at node 40 (Figure 2.1). 
This location would result in the weakest signal and necessarily correspond to the most difficult 


NAS A/TM— 2004-2 13118 


20 



location for FOD event detection. At 3.44 seconds, a 0.5 lbm foreign object hits the fan disk with a 
speed of 300 mph (parallel to the Z-axis in Figure 2.3) at a radius of 20 inches. The event is 
modeled as a pulse of width 0.04 seconds, with a magnitude determined using the method 
presented in Section 2.2. A fan disk eccentricity of 0.001 inches is assumed. As shown in 
Figure 4.1a , the event is barely noticeable in the time trace for a noise-free situation. However 
accelerometer signals on in-service engines are typically noisy. In addition to “process noise” i.e., 
random motion of the aircraft due to wind gusts and compensating maneuvers being transmitted 
from the airframe to the engine, there is a significant amount of sensor noise which may mask the 
occurrence of a FOD event. Figure 4.1b shows the accelerometer signal at location 40, with a 
signal-to-noise ratio of approximately 3.5 (11 dB). Observation of the signal shows no well-defined 
point in time at which one would identify a FOD event occurrence. 



Time (sec) 


50 



0 I 1 1 1 I L_ 

3.3 3.35 3.4 3.45 3.5 3.55 

Time (sec) 


Figure 4. 1 . Bearing 40 accelerometer output signal corresponding to a FOD event, 
a) Noise free system response, b) System response with process and sensor noise. 
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4.1 FOD event Detection using Fourier Analysis of Rotor Bearing Vibration Signals 


Vibration signal analysis is an invaluable tool for identifying mechanical component problems. 
The most utilized tool in vibration analysis i.e., the Fourier transform, essentially assumes that the 
signals being analyzed are of infinite duration [18,19]. When the continuous-time Fourier transform 


Y(F)= P y(t)e -j2;rFt dt 

J — oo 


4.1 


or, more practically, its discrete-time counterpart 


n= oo 

Y(f)= X y(n)e“ j2;rfn 


n=-oo 


4.2 


is applied to a signal of finite duration, spectral leakage effects can significantly reduce the 
resolution of the corresponding power (frequency) spectrum given by 


|Y(fJ 


2 


I y(n)e -j2jrfn 


2 


n= -oo 


4.3 


Windowing of the original signal aids in mitigating the effects of finite signal length (spectral 
leakage), however if the component of the signal being analyzed has “compact support” (e.g., an 
impulse due to a FOD-event induced structural response), windowing of the signal may not 
significantly aid in identifying the occurrence of the event in time due to the Fourier transform 
basis functions i.e., sinusoids, being themselves functions of infinite support and limited in their 
ability to decompose such a finite-duration signal. The short-time Fourier transform (and the 
corresponding time-dependant power spectrum) uses a windowed version of y(n) to compute Y(f) 
in an attempt to localize a change in a signal over time [19]. The time-dependant power spectrum 
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for the noise-free signal (Figure 4.2 a) shows a significant change at approximately 3.3 seconds, 
which disappears as time progresses. The time-dependant power spectrum of the noisy signal 
(Figure 4.2 b) shows no noticeable change in time, with the effect of the FOD event “buried” in the 
spectra due to the noise. Thus, there is a need for an alternative signal analysis technique for FOD 
event detection. 



Figure 4.2. Power spectra of bearing 40 accelerometer output signal corresponding to a 
FOD event, a) Noise free system response, b) System response with process and sensor 
noise. 
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4.2 FOD event Detection using Wavelet Analysis of Rotor Bearing Vibration Signals 


Over the past decade the discrete time wavelet transform (DTWT) has been applied to a wide 
range of signal analysis problems e.g., de-noising of signals as well as time localization and 
reconstruction of short duration changes [8]. Because of these characteristics, the DTWT is 
considered to be a viable candidate for identifying a short duration change in a bearing vibration 
signal due to a FOD event corrupted by noise. 

The DTWT is the discrete-time counterpart of the continuous-time wavelet transform 

CWT(a, r) = J y(t) jdt 4.4 

where v P(t) is the basic or mother wavelet and v P((t - x)/a)/Va are the wavelet basis functions which 
are scaled via the parameter a (by compressing or stretching the mother wavelet) and shifted in 
time by x . The operation performed in equation 4.4 on the original signal y(t) may be interpreted 

as the cross correlation of the signal y(t) with v F(t/a)/Va , shifted by x/a [18]. Thus, the CWT 

computes the component of y(t) that is similar to v P(t/a)/Va . If little (or no) similarity exists 
between the two, then the CWT will be small (or zero). Larger values of CWT indicate better 
correlation. Figure 4.3 presents an overview of the Wavelet Transform process with multiple scale 
decomposition. The wavelet family used is the Daubechies family. 

Typically, the resulting data is presented in terms of time (or k if in discrete time) and scale 
(i.e., coefficient) of the chosen wavelet transform. “Pseudo” frequency components of y(t), 
corresponding to the scale of the constituent wavelets, F m ,are obtained by the following 
expression [20] 



m 


4.5 
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where m is the scale, T is the sampling interval, and F c is the wavelet center frequency calculated 

by taking the maximum of the modulus of the wavelet’s Fourier transform. 

Most of the functions for 'F(t) have been developed over the past two decades with the Haar 

wavelet being the first documented function to be used in this context[19]. Choice of 'F(t) is 

heavily dependant on the nature of the signal (or a localized feature of a signal) being analyzed. For 

example, if a vibration signal were to be analyzed to determine whether a FOD event had occurred, 

one would expect to see an impulse-like localized feature embedded in a normally noisy sinusoidal 

signal. Intuitively, one would choose a wavelet that had an impulse-like character as opposed to 

one that resembled a step function, as is the case with the Haar wavelet mentioned previously. 

Inverse Wavelet transform 

gives signal decomposition in terms 

of time and scale (pseudo frequency) 



Figure 4.3 Example bearing vibration signature and corresponding wavelet 
decomposition accompanying a FOD event 

In discrete-time the wavelet transform takes the form 

m 

2 

DTWTD(m,n) = a 0 ^y(k)y(a 0 -m k - nx 0 ) 4.6 

k 
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where m is the scale (degree of dilation) and n corresponds to the net translation in time (in terms 
of the number of time x 0 ) of the wavelet at a specific scale. For scale m=0, the coefficients of the 
so-called scaling function, from which all of the wavelets in a given family are related, are 


DTWTA(0, n) = £ y(k)cp(k - nx 0 ) 4.7 

k 

tp(k - nx 0 )is typically determined via iteration using the following dilation equation 


N-l 

(p(k) = ^c n (p(2k-nx 0 ) 


n=0 


4.8 


where N is the number of terms in (p(k) . The scaling function cp(k) determines the character of 
the wavelet. The characteristic of interest in the signal being analyzed typically dictates the form of 
the wavelets used. A wavelet \|/(k) is related to it’s associated cp(k) via 


N-l 


V(k) = Yj (' 1 )' 1 c „<P(2k + nx 0 - N + l) 


n=0 


4.9 


Equation 4.9 is also a dilation equation with \|/(k) referred to as a dilation wavelet [19]. The 
scaling function, it’s dilates, the corresponding wavelets and their dilates are orthogonal to one 
another and therefore constitute a set of functions by which an arbitrary function may be built [19], 
completely analogous to the Fourier transform and series. The base value of scaling parameter, a 0 , 

and the time shift, x 0 , are typically set equal to 2 and 1 respectively for computational efficiency, 

analogous to the Fast Fourier Transform (FFT). 

The convolution of Equations 4.6 and 4.7 is commutative i.e., 
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4.10 


DTWTA(0,n) = - nx G ) 

k 

m 

DTWTD(m, n) = a 0 " y X Y(k)y(a 0 “ m k - nx 0 ) 4.11 

k 

Thus, the discrete-time wavelet transform may be viewed as a filtering operation with 
cp(k) containing the coefficients of a low -pass Finite Impulse Response (FIR) filter and \|/(k) the 
coefficients of a high-pass FIR filter. The coefficients corresponding to the low-frequency portion 
of the input signal i.e., approximation coefficients (Equation 4.10), and the high-frequency detail 
coefficients (Equation 4.1 1) at a given scale in a wavelet decomposition are provided by successive 
low pass and high pass filtering operations respectively. Reconstruction of the approximations and 
details at a specified scale m is performed via the corresponding inverse wavelet transform 


D m (k) = X DTWTD(m, n)y m>n (k) 

n 

4.12 

A m (k) = X DTWTA(m, n) cp m n (k) 

4.13 


n 


and is implemented via banks of reconstruction or synthesis filters. Figure 4.4 presents a 2-scale 
DTWT-based signal analysis and synthesis using an FIR quadrature mirror filter bank. Locations in 
the filter bank correspond to the signal characteristics of interest (details or approximations). 
Scaling is accomplished by downsampling after high or low pass filtering. Successive upsampling 
and filtering performs the reconstruction. When the reconstruction filter coefficients correspond to 
the quadrature mirror filter of the decomposition filter, perfect signal reconstruction is guaranteed 
[21]. This is the practical implementation of the discrete-time wavelet transform algorithm 
developed by Mallat [22, 23, 24]. The information gained from the decomposition depends on the 
application e.g., a signal may be “denoised” by omitting the detail from one or more scales during 
the reconstruction. As well, abrupt changes in the signal e.g., a FOD event-induced accelerometer 
response, might be detected by observing one or more of the approximations. 
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Figure 4.5 presents a wavelet analysis of the signal analyzed in Section 4.1, the noise-free 
bearing accelerometer signal at location 40. The Daubechies 8 wavelet was chosen (a subjective 
choice) for the analysis due to its correspondence to the signal characteristic of interest. Observing 
the approximation at scale 8 (Figure 4.5b) demonstrates the wavelet’s ability to determine the 
location in time of the event. Several other approximations and details at various scales (not shown) 
also revealed the occurrence of a FOD event. The analysis of the noisy signal at the same scale’s 
approximation (Figure 4.6), however, illustrates that the wavelet chosen (as well as several other 
wavelets tested, not presented here), may have difficulty highlighting a subtle change buried in a 
low signal-to-low noise ratio (SNR = 3.5) signal. Additional “conditioning” of the decomposed 
signal is required. 



Scale 2 Approximation 



Low Pass Reconstruction FIR 
Filter 



Low Pass Decomposition FIR 
Filter 

High Pass Reconstruction FIR 
Filter 



High Pass Decomposition FIR 
Filter 



Downsample 
by a Factor of 2 



Upsample 
by a Factor of 2 


Figure 4.4. Example of a 2-scale DTWT-based signal analysis and synthesis using an FIR 
quadrature mirror filter bank. 
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Figure 4.5. Wavelet analysis of noise-free bearing accelerometer signal at location 40. 
a) accelerometer signal, b) Corresponding Daubechies 8 wavelet inverse transform, scale 8 
approximation. 




a) accelerometer signal, b) Corresponding Daubechies 8 wavelet inverse transform, scale 8 
approximation. 
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Abrupt changes in an otherwise smooth signal will correspond to a localized change in the 
signal derivative with respect to time (or some space-related dependant variable depending on the 
application). The more abrupt the change, the larger the local value of the derivative with respect to 
time. Additionally, for the same change the value of higher order derivatives 


d n y(t) 

dt n 


4.14 


will provide a frequency-dependant gain i.e., the more abrupt the signal change, the larger the gain. 
This quality of the derivative is more easily observed in the frequency domain by taking the 
Laplace transform of Equation 4.14 and setting s = jco 


T kV) 

{ dt n 


s=jco 


(jo) n y(jco) 


4.15 


The corresponding frequency-dependant magnitude is 

MAG^( jco) n y(jco)j = co n MAG(y( jco)) 


4.16 


Thus, multiple differentiation in the time domain provides an “s 11 ” multiplication factor in the 
frequency domain. The higher frequency characteristics of a signal (e.g., a FOD-event would be 
classified as a high-frequency event) would have a significantly higher gain relative to the low- 
frequency components. It would seem that in a noise-free signal, the occurrence of a FOD-event 
resulting in an impulsive input to the turbofan rotor of the present study could easily be identified 
using this technique. Multiple differentiation of noisy signals is typically avoided in most signal 
processing applications due to the amplification provided by Equation 4.16. Indeed, multiple 
differentiation of the signal shown in Figure 4.6a would make identification of a FOD-event 
virtually impossible, primarily due to amplification of the noise. 


NAS A/TM— 2004-2 13118 


30 



Mallat [25] proposed a technique for identifying a wavelet specifically designed for edge 
detection in computer vision systems. A signal (or 2-D image in the original application) is passed 
through a smoothing function (e.g., a low pass filter) and differentiated multiple times. The 
combination of smoothing and differentiation results in a wavelet customized to the application. 
For this research, the technique was adapted to identify the FOD event induced short-time change 
in accelerometer signals. The smoothing filter chosen for the present investigation is the 
Daubechies-8 wavelet decomposition. As shown in Figure 4.7, the accelerometer signal is passed 
through an approximation filter bank three times and subsequently passed through a detail filter. 
This provides 




Low Pass Reconstruction 
FIR Filter 


Lj} Low Pass Decomposition 

1 FIR Filter 

High Pass Reconstruction 
1 FIR Filter 



High Pass Decomposition 
FIR Filter 



Downsample 
by a Factor of 2 



Upsample 
a Factor of 2 


Figure 4.7. Wavelet transform-based vibration feature extraction implemented using 
an analysis and synthesis FIR filter bank. 
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a means to focus on a subband of the original signal where the feature of interest is thought to lie, 
the net effect being a low-pass filtering or smoothing operation. The result is differentiated three 
times to produce the wavelet-based feature of interest. Figure 4.8 presents the wavelet feature 
extracted from the accelerometer signal shown in Figure 4.6, using the custom wavelet described 
above. The feature extracted is the spike shown in Figure 4.6b. The exact time of the event is 
precisely determined using this technique. Indeed, the event is modeled by imposing a pulse, whose 
magnitude is a function the of foreign object impact characteristics mentioned previously, on the 
fan disk. The two spikes shown in Figure 4.8 result from the rising and descending parts of that 
pulse “filtered” through the mechanical (i.e., rotor-bearing) system. The event input profile to the 
model (i.e., the pulse input described above) was considered to be appropriate for a significant, 
structurally damaging event i.e., where the foreign object absorbs virtually none of the energy of 
impact, resulting in complete and immediate transfer of energy to the fan disk. This would, for 
example, be the case for hard objects such as ice. For other scenarios, pulses of lower magnitude or 
with non-infinite rising (or descending) slopes may require higher-order differentiation for 
detection. 

In summary, the result of the following sequence of operations: 

• wavelet decomposition of a noisy ( e.g., bearing accelerometer) signal 

• reconstruction of selected components (approximations or details at a selected scales) of the 
signal 

• multiple differentiation of the reconstructed components 

can precisely time-locate an abrupt, high-frequency, change in the signal. However, choice of too 
high a scale may result in the high-frequency event being filtered out of the signal. 
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Figure 4.8. Wavelet feature extracted from noisy bearing accelerometer signal at location 
40. a) accelerometer signal, b) Corresponding conditioned Daubechies 8 wavelet inverse 
transform, scale 8 approximation. 


5.0 Conclusion 


This paper describes the development of a wavelet-based feature (WF) specifically designed to 
identify the occurrence of FOD-events in gas turbine engines. The WF is extracted from rotor 
bearing accelerometer signals and is shown to be robust in the presence of significant noise. The 
technique is intended ultimately be combined with analytical measurements (e.g., Kalman filter 
thermal/health parameter estimates) for FOD-event detection via information fusion from these 
(and perhaps other) sources. It is envisioned that fusion of structural and thermal performance 
information will provide a pilot, or an autonomous vehicle, more confidence in a FOD-event 
diagnosis thus, enabling a corrective action to be based on a more informed decision. 
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Due to the lack of high-frequency FOD-event test data in the open literature, a reduced-order 
turbofan structural model (ROM) was synthesized from a finite element model to support the 
investigation. The wavelet-based FOD-event detection scheme developed is “data driven” (i.e., not 
model-based), thus precise correspondence of the ROM to the FEM (which itself would probably 
not represent an actual engine with complete fidelity) is not essential. The ROM is considered to be 
“phenomenologically accurate,” and shown to provide acceptable fidelity with regards to the 
application i.e., real-time diagnostic system development. 

Using the state-space ROM, a sensor placement analysis was performed to determine the 
optimal accelerometer location for FOD-event detection. To demonstrate the robustness of the 
technique, testing was performed using the least optimal bearing location. For comparison, analysis 
of noise-free and noisy accelerometer signals was performed using a standard wavelet 
decomposition and shown to be inadequate for detecting FOD events. In the presence of significant 
noise (SNR=3.5, 1 1 dB), precise location of the FOD event in time was obtained using the custom 
wavelet feature developed. Choice of the smoothing wavelet used ( Daubechies 8) was subjective, 
and based on knowledge of the signal characteristic of interest i.e., a FOD event-induced abrupt 
change in an otherwise noisy signal. 

Future work should include determining the optimal wavelet (and corresponding scale) to apply 
based on current information theoretic techniques such as maximum entropy. Additionally, 
combination of structural and thermal performance information could provide a robust diagnosis-of 
foreign object damage in specific engine components i.e., extending this research to the Fault 
Detection Isolation (FDI) problem. Continued development and verification of the feature 
presented would require high-frequency accelerometer data (i.e., sensor bandwidth greater than 10 
kHz) from FOD event engine tests at multiple operating conditions. 

6.0 Nomenclature 

A m Reconstructed approximation component of wavelet-analyzed signal at scale m 

M Mass Matrix 

C Damping Coefficient Matrix 

D m Reconstructed detail component of wavelet-analyzed signal at scale m 


NAS A/TM— 2004-2 13118 


34 



DTWTD 

DTWTA 

K 

G 

Cxy>K X y 

a,P 

(0 

F,f 

F 

^0 > U 
m 

n 

T|(t) 

Y(f) 

Vm,n 


q(t) 

<P\ m,n 


Discrete-time wavelet transform detail coefficient 

Discrete-time wavelet transform approximation coefficient 

Stiffness, Stiffness Matrix 
Gyroscopic Matrix 

Quantities with coupling between, e.g., the x and y degrees of freedom 
Proportional damping coefficients 
angular frequency (rad/sec) 

Frequency (Hz) 

Force 

Moments of inertia about the X and Y axes respectively 
Dyadic scale of wavelet-decomposed signal 

Wavelet time delay in number of samples 

Displacement in modal coordinates 
Discrete Fourier Transform of y(n) 

Wavelet at scale m and translation n 
Displacement in physical coordinates 

Wavelet scaling function at scale m and translation n 
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